1 Introduction

The objective of this notebook is to perform a basic quality control of the mapping performed with Cellranger. Cell Ranger ARC analyzes data generated by the Chromium Single Cell Multiome ATAC + Gene Expression. Should note here, we will talk about features that can refer to a gene (declared from the hg38 reference transcriptome) or a peak (determined by the Peak Calling algorithm provided by Cell Ranger).

The expected values (EV) included in the notebook are obtained from this technical_note

2 Pre-processing

2.1 Load packages

library(ggplot2)
library(ggpubr)
library(plyr)
library(reshape2)
library(data.table)
library(knitr)
library(kableExtra)
library(ggrepel)
library(plotly)

2.2 Parameters

path_to_multiome <- here::here("multiome/results/tables/cellranger_mapping/cellranger_mapping_metrics_multiome.csv")
path_to_proj_metadata <- here::here("multiome/1-cellranger_mapping/data/tonsil_atlas_metadata_multiome.csv")

2.3 Functions

barplot_data <- function(multiome_melted){
  multiome_melted$value <- round(multiome_melted$value, 2)
  ggbarplot(multiome_melted,
  x = "library_name",
  y = "value",
  fill = "gray70",
  ggtheme = theme_pubr(x.text.angle = 45,legend = "none")) + 
    facet_wrap(ncol = 1,scales = "free_y", ~variable) 
}

correlation_plot <- function(multiome_df, variable1, variable2)
ggscatter(qc_samples, x = variable1, y = variable2,
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   ) + stat_cor(method = "pearson")

2.4 Load data

qc_samples <- read.csv(path_to_multiome)
metadata <- read.csv(path_to_proj_metadata)
qc_samples$library_name <- revalue(
  qc_samples$Sample.ID,
  c(
    "ulx1v6sz_8a2nvf1c" = "BCLL_8_T_1",
    "wdp0p728_jf6w68km" = "BCLL_8_T_2",
    "co7dzuup_xuczw9vc" = "BCLL_9_T_1",
    "qmzb59ew_t11l8dzm" = "BCLL_9_T_2",
    "pd9avu0k_kf9ft6kk" = "BCLL_14_T_1",
    "vuuqir4h_wfkyb5v8" = "BCLL_14_T_2",
    "admae8w2_89i88tvv" = "BCLL_15_T_1",
    "sr20954q_yiuuoxng" = "BCLL_15_T_2",
    "kmbfo1ab_ie02b4ny" = "BCLL_2_T_1",
    "ryh4el3i_biv0w7ca" = "BCLL_2_T_2",
    "bs2e7lr7_mdfwypvz" = "BCLL_2_T_3"
    
  )
)
knitr::kable(qc_samples, row.names = TRUE) %>%
  kable_styling(
    bootstrap_options = c("striped"),
    full_width = T,
    font_size = 12
  ) %>%
  scroll_box(height = "300px")
Sample.ID Pipeline.version Genome Estimated.number.of.cells ATAC.Confidently.mapped.read.pairs ATAC.Fraction.of.genome.in.peaks ATAC.Fraction.of.high.quality.fragments.in.cells ATAC.Fraction.of.high.quality.fragments.overlapping.TSS ATAC.Fraction.of.high.quality.fragments.overlapping.peaks ATAC.Fraction.of.transposition.events.in.peaks.in.cells ATAC.Mean.raw.read.pairs.per.cell ATAC.Median.high.quality.fragments.per.cell ATAC.Non.nuclear.read.pairs ATAC.Number.of.peaks ATAC.Percent.duplicates ATAC.Q30.bases.in.barcode ATAC.Q30.bases.in.read.1 ATAC.Q30.bases.in.read.2 ATAC.Q30.bases.in.sample.index.i1 ATAC.Sequenced.read.pairs ATAC.TSS.enrichment.score ATAC.Unmapped.read.pairs ATAC.Valid.barcodes Feature.linkages.detected GEX.Fraction.of.transcriptomic.reads.in.cells GEX.Mean.raw.reads.per.cell GEX.Median.UMI.counts.per.cell GEX.Median.genes.per.cell GEX.Percent.duplicates GEX.Q30.bases.in.UMI GEX.Q30.bases.in.barcode GEX.Q30.bases.in.read.2 GEX.Reads.mapped.antisense.to.gene GEX.Reads.mapped.confidently.to.exonic.regions GEX.Reads.mapped.confidently.to.genome GEX.Reads.mapped.confidently.to.intergenic.regions GEX.Reads.mapped.confidently.to.intronic.regions GEX.Reads.mapped.confidently.to.transcriptome GEX.Reads.mapped.to.genome GEX.Reads.with.TSO GEX.Sequenced.read.pairs GEX.Total.genes.detected GEX.Valid.UMIs GEX.Valid.barcodes Linked.genes Linked.peaks library_name
1 qmzb59ew_t11l8dzm cellranger-arc-1.0.0 GRCh38 5530 0.90964 0.05144 0.90915 0.44015 0.77105 0.75653 57807.81 8271 0.02425 132440 0.70263 0.91053 0.96303 0.96083 NA 319677205 11.26773 0.01051 0.97903 151554 0.90220 39397.02 4847.0 2142 0.68264 0.96989 0.97308 0.95016 0.22107 0.37121 0.91608 0.05495 0.48992 0.63371 0.96232 0.05575 217865541 28636 0.99949 0.94706 11028 45856 BCLL_9_T_2
2 co7dzuup_xuczw9vc cellranger-arc-1.0.0 GRCh38 5572 0.90972 0.04413 0.85491 0.42580 0.74570 0.72928 55415.02 8776 0.02305 115988 0.66436 0.89602 0.95478 0.95245 NA 308772502 10.55462 0.00945 0.97738 147451 0.91009 35104.20 4322.0 1997 0.68069 0.97089 0.97399 0.94723 0.22482 0.36198 0.92291 0.05291 0.50802 0.63881 0.96428 0.04315 195600575 28243 0.99947 0.94793 10657 42827 BCLL_9_T_1
3 wdp0p728_jf6w68km cellranger-arc-1.0.0 GRCh38 7861 0.90808 0.03255 0.76271 0.40009 0.68396 0.66529 31841.03 6958 0.00803 98142 0.49401 0.91650 0.96424 0.96264 NA 250302362 9.61085 0.00692 0.98580 157487 0.92367 24671.96 4557.0 2117 0.52518 0.97029 0.97346 0.95037 0.22588 0.38523 0.92960 0.05497 0.48940 0.64239 0.96983 0.05795 193946270 28847 0.99949 0.94805 11387 38686 BCLL_8_T_2
4 ulx1v6sz_8a2nvf1c cellranger-arc-1.0.0 GRCh38 7181 0.90888 0.03716 0.80033 0.42550 0.72507 0.70834 31305.98 7110 0.00855 108000 0.49722 0.91854 0.96808 0.96549 NA 224808216 10.26944 0.00778 0.98630 170819 0.92461 30537.81 5236.0 2318 0.55138 0.96930 0.97212 0.95113 0.22932 0.37099 0.92631 0.05645 0.49887 0.63400 0.97059 0.05456 219291981 29095 0.99950 0.94705 12085 42358 BCLL_8_T_1
5 admae8w2_89i88tvv cellranger-arc-1.0.0 GRCh38 5788 0.89701 0.02220 0.60932 0.45470 0.64987 0.62956 39558.28 7394 0.00290 74187 0.42259 0.90235 0.95770 0.95480 NA 228963303 10.45955 0.01080 0.98068 139918 0.87102 24146.67 3476.5 1713 0.56009 0.97060 0.97394 0.95230 0.17797 0.43160 0.91885 0.06388 0.42336 0.67030 0.97442 0.06432 139760911 27740 0.99953 0.95331 8430 30643 BCLL_15_T_1
6 sr20954q_yiuuoxng cellranger-arc-1.0.0 GRCh38 5845 0.89716 0.02227 0.61257 0.46364 0.66085 0.64027 39362.91 7073 0.00300 73885 0.42324 0.90302 0.95704 0.95494 NA 230076230 10.55038 0.01108 0.98122 140366 0.88381 30731.03 4005.0 1895 0.60792 0.96944 0.97246 0.95280 0.17746 0.43255 0.91632 0.06370 0.42007 0.66831 0.97285 0.06891 179622850 28285 0.99953 0.95371 8855 30571 BCLL_15_T_2
7 bs2e7lr7_mdfwypvz cellranger-arc-1.0.0 GRCh38 7460 0.88155 0.01497 0.44372 0.41225 0.52801 0.50781 26579.67 4298 0.01152 54687 0.26634 0.89865 0.95679 0.95465 NA 198284372 9.51567 0.00723 0.97758 19393 0.90605 21134.62 3354.0 1586 0.55768 0.96849 0.97215 0.94678 0.12086 0.35940 0.91054 0.07914 0.47199 0.70309 0.97493 0.11828 157664248 28650 0.99964 0.96253 5099 10548 BCLL_2_T_3
8 pd9avu0k_kf9ft6kk cellranger-arc-1.0.0 GRCh38 6919 0.91678 0.05087 0.91946 0.44305 0.73416 0.71805 24739.50 7714 0.00444 123616 0.44733 0.88822 0.94466 0.94412 NA 171172604 10.70481 0.01136 0.98087 144699 0.91809 28372.47 3828.0 1736 0.67548 0.96954 0.97278 0.95094 0.17905 0.41910 0.91917 0.06488 0.43519 0.66884 0.96896 0.06737 196309095 28051 0.99954 0.95416 9922 45191 BCLL_14_T_1
9 vuuqir4h_wfkyb5v8 cellranger-arc-1.0.0 GRCh38 6103 0.91469 0.04771 0.90557 0.43784 0.72039 0.70473 25950.80 7821 0.00411 120001 0.42968 0.89776 0.95826 0.95511 NA 158377744 10.74407 0.01220 0.97978 123546 0.90995 23631.59 3561.0 1659 0.62458 0.97037 0.97379 0.94388 0.17702 0.41979 0.92120 0.06485 0.43656 0.67284 0.96767 0.08440 144223620 27459 0.99958 0.95461 8947 40870 BCLL_14_T_2
10 kmbfo1ab_ie02b4ny cellranger-arc-1.0.0 GRCh38 10837 0.87842 0.01395 0.46442 0.36849 0.45577 0.43771 23010.08 3623 0.00997 51351 0.25441 0.90147 0.95567 0.95404 NA 249360235 8.97706 0.00718 0.98139 15409 0.92051 14295.57 2425.0 1291 0.49744 0.96769 0.97084 0.94821 0.12348 0.35548 0.90051 0.07662 0.46840 0.69285 0.97212 0.11530 154921110 29143 0.99959 0.96211 4788 9182 BCLL_2_T_1
11 ryh4el3i_biv0w7ca cellranger-arc-1.0.0 GRCh38 7910 0.88602 0.01503 0.45653 0.40099 0.51282 0.49108 17759.52 3088 0.01001 53499 0.21629 0.88597 0.94668 0.94405 NA 140477823 9.16290 0.00793 0.97915 13469 0.90946 20475.87 3298.5 1584 0.54781 0.97118 0.97343 0.94894 0.11372 0.35705 0.89950 0.07922 0.46323 0.69910 0.96995 0.11509 161964102 28869 0.99957 0.96247 4640 8671 BCLL_2_T_2

3 Multiome Metrics Definitions

3.1 Estimated Number of Cells

Let us start by plotting the estimated number of cells per library. Note this this number will differ a lot from the final number of cells after applying future QC filters.

barplot_data(melt(qc_samples[,c("library_name","Estimated.number.of.cells")])) + coord_flip()

print(paste("Estimated number of cells for all the samples:", sum(qc_samples$Estimated.number.of.cells)))
## [1] "Estimated number of cells for all the samples: 77006"

3.2 Joint Metrics

Here, we can visualize the number of the features (peaks and genes) and feature linkages detected. Note that the feature linkage is calculated taking in account genes with its proximal peaks and between pairs of proximal peaks across the genome (based on correlation).

We can see how three samples BCL_2_T_3, BCL_2_T_1, BCL_2_T_2 have a low number of detected features, maybe caused by a low number of recovered nuclei These samples could be classified a priori as low quality samples.

barplot_data(melt(qc_samples[,c("library_name",
                                "Linked.genes",
                                "Linked.peaks")]))

linkage_melt <- melt(qc_samples[,c("library_name",
                                   "Feature.linkages.detected")])
barplot_data(multiome_melted = linkage_melt)

4 Chromatin Accessibility Metrics Definitions

4.1 Sequencing

Our first objective is to evaluated the quality and the quantity of our sequenced libraries prior to mapping.

4.1.1 Sequenced read pairs

Here, you can see the total number of sequenced read pairs assigned to the Multiome ATAC library. The optimal number of read pairs are 25,000 per cell taking in account the technical note.

qc_samples$ATAC.Mean.raw.read.pairs.per.cell <- round(qc_samples$ATAC.Mean.raw.read.pairs.per.cell,2)
ggbarplot(qc_samples,
  x = "library_name",
  y =  "ATAC.Mean.raw.read.pairs.per.cell",
  label = TRUE,
  fill = "gray70",
  ggtheme = theme_pubr(x.text.angle = 45,legend = "none")) + 
  geom_hline(yintercept  = 25000, linetype='dotted', col = 'black')

Here, you can see the quality (leveraging the “Q30” variables) of sequenced read pairs assigned to the Multiome ATAC library.

multiome_melted <- melt(qc_samples[,c("library_name", "ATAC.Q30.bases.in.barcode",
"ATAC.Q30.bases.in.read.1",        
"ATAC.Q30.bases.in.read.2")])
          
ggplot(multiome_melted, aes(variable, value)) +
  geom_boxplot(fill = NA, outlier.shape = NA) +
  geom_jitter(width = 0.1, height = 0) +
  geom_text_repel(data = multiome_melted, aes(label = library_name, fill = "black")) +
  theme_minimal() +
  theme(legend.position = "none",
      axis.title = element_text(size = 10),
      axis.text = element_text(size = 10, color = "black"))

4.1.2 Valid barcodes

Fraction of read pairs with barcodes that match the whitelist after error correction. A value higher than 85% represent a high quality library.

ggbarplot(qc_samples,
  x = "library_name",
  y =  "ATAC.Valid.barcodes",
  label = TRUE,
  fill = "gray70",
  ggtheme = theme_pubr(x.text.angle = 45,legend = "none")) + 
  geom_hline(yintercept  = 0.85, linetype='dotted', col = 'black')

4.1.3 Percent of duplicates

The percentatge of the duplicated data is correlated with the library size and the sequencing depth.

correlation_plot(qc_samples,"ATAC.Percent.duplicates","ATAC.Sequenced.read.pairs")

correlation_plot(qc_samples,"ATAC.Percent.duplicates","ATAC.Number.of.peaks")

4.2 Cells

4.2.1 Median number of fragments per cell

High-Quality Fragment: A read-pair with mapping quality > 30, that is not chimerically mapped, has a valid 10x barcode, and maps to any nuclear contig (not mitochondria) that contains at least one gene.

qc_atac_read_pairs <- melt(qc_samples[,c("library_name","ATAC.Mean.raw.read.pairs.per.cell",                      "ATAC.Median.high.quality.fragments.per.cell")])
barplot_data(qc_atac_read_pairs)

qc_atac_fraction <- melt(qc_samples[,c("library_name","ATAC.Fraction.of.high.quality.fragments.in.cells","ATAC.Fraction.of.transposition.events.in.peaks.in.cells")])
barplot_data(qc_atac_fraction)

4.2.2 Targeting Metrics

The targetting metrics can be summarized by these 4 main score:

4.2.2.1 Total number of peaks on primary contigs.

This number presents a high correlation with the sequencing depth specifically with the high quality fragments.

ggbarplot(qc_samples, "library_name", "ATAC.Number.of.peaks",
  label = TRUE,
  fill = "gray70",
  ggtheme = theme_pubr(x.text.angle = 45,legend = "none"),
  position = position_dodge(0.9))

 correlation_plot(qc_samples,"ATAC.Number.of.peaks",
 "ATAC.Median.high.quality.fragments.per.cell")

4.2.2.2 Fraction of high quality fragments that overlapp in peaks.

The fraction of high quality fragments in cell are expected to be higher 40%. The fraction of transposition events that fall within peaks > 25%

multiome_melted <- melt(qc_samples[,c("library_name", "ATAC.Fraction.of.high.quality.fragments.in.cells",
"ATAC.Fraction.of.high.quality.fragments.overlapping.peaks",        "ATAC.Fraction.of.high.quality.fragments.overlapping.TSS")])
                 
multiome_melted$value = round(multiome_melted$value,2)                   
ggbarplot(multiome_melted, "library_name", "value",
  fill = "variable", color = "variable", palette = "Paired",
  label = TRUE,
  ggtheme = theme_pubr(x.text.angle = 45),
  position = position_dodge(0.9))

4.2.2.3 Fraction of genome in peaks

This fraction is quite low in all samples. We do not have a reference value to be able to compare them. However, the examples we reviewed had a value of around 2%.

ggplot(melt(qc_samples[,c("library_name",
                          "ATAC.Fraction.of.genome.in.peaks")]), aes(variable, value)) +
  geom_boxplot(fill = NA, outlier.shape = NA) +
  geom_jitter(width = 0.1, height = 0) +
  theme_minimal() +
  theme(legend.position = "none",
      axis.title = element_text(size = 10),
      axis.text = element_text(size = 10, color = "black"))

4.2.2.4 Transcription Start Site (TSS)

It is expected to obtain a large enrichment around TSS (> 5% EV, red dashed line), as these regions are known to show a high degree of chromatin accessibility compared to the flanking regions.

data_tss <- qc_samples[, c("library_name", "ATAC.TSS.enrichment.score")]
data_tss.melt <- melt(data_tss)
barplot_data(data_tss.melt) + geom_hline(yintercept = 5, linetype = 2, color = "red")

5 Gene Expression Metrics Definitions

5.1 Median genes per cell

median_genes_gg <- melt(qc_samples[,c("library_name","GEX.Median.genes.per.cell")])
barplot_data(median_genes_gg)

5.1.1 Sequenced read pairs

Here, you can see the total number of sequenced read pairs assigned to the Multiome Gene Expression library. The optimal number of read pairs are ~20,000 per cell taking in account the technical note.

qc_samples$GEX.Mean.raw.reads.per.cell <- round(qc_samples$GEX.Mean.raw.reads.per.cell,2)
ggbarplot(qc_samples,
  x = "library_name",
  y =  "GEX.Mean.raw.reads.per.cell",
  label = TRUE,
  fill = "gray70",
  ggtheme = theme_pubr(x.text.angle = 45)) + 
  geom_hline(yintercept  = 20000, linetype='dotted', col = 'black')

Here, you can see the quality (leveraging the “Q30” variables) of sequenced read pairs assigned to the Multiome Gene Expresion library.

multiome_melted <- melt(qc_samples[,c("library_name", "GEX.Q30.bases.in.barcode",
"GEX.Q30.bases.in.UMI",        
"GEX.Q30.bases.in.read.2")])

ggplot(multiome_melted, aes(variable, value)) +
  geom_boxplot(fill = NA, outlier.shape = NA) +
  geom_jitter(width = 0.1, height = 0) +
  geom_text_repel(data = multiome_melted, aes(label = library_name, fill = "black")) +
  theme_minimal() +
  theme(legend.position = "none",
      axis.title = element_text(size = 10),
      axis.text = element_text(size = 10, color = "black"))

5.1.2 Valid barcodes

Fraction of read pairs with barcodes that match the whitelist after error correction. A value higher than 80% represent a high quality library,

ggbarplot(qc_samples,
  x = "library_name",
  y =  "GEX.Valid.barcodes",
  label = TRUE,
  fill = "gray70",
  ggtheme = theme_pubr(x.text.angle = 45,legend = "none")) + 
  geom_hline(yintercept  = 0.85, linetype='dotted', col = 'black')

5.1.3 Percent of duplicates

The percentatge of the duplicated data is correlated with the library size and the sequencing depth.

correlation_plot(qc_samples,"GEX.Percent.duplicates","GEX.Sequenced.read.pairs")

correlation_plot(qc_samples,"GEX.Percent.duplicates","GEX.Total.genes.detected")

5.2 Mapping

Secondly, we will assess the quality of cellranger’s mapping by comparing the percentage of reads mapping to the genome, intergenic regions, intronic and exonic regions across libraries. Reads mapping to intergenic regions suggest contamination of ambient DNA, while reads mapping to intronic regions may come from pre-mRNAs or mature splice isoforms that retain the introns.

The fraction of sequenced reads that map to a unique gene in the transcriptome is expected to be higher than the 50%.

qc_samples$GEX.Fraction.of.transcriptomic.reads.in.cells <- round(qc_samples$GEX.Fraction.of.transcriptomic.reads.in.cells,2)
ggbarplot(qc_samples,
  x = "library_name",
  y =  "GEX.Fraction.of.transcriptomic.reads.in.cells",
  fill = "gray70",
  label = TRUE,
  ggtheme = theme_pubr(x.text.angle = 45)) + 
  geom_hline(yintercept  = 0.50, linetype='dotted', col = 'black')

multiome_melted <- melt(qc_samples[,c("library_name", "GEX.Reads.mapped.to.genome",  "GEX.Reads.mapped.confidently.to.genome")])
                
multiome_melted$value <- round(multiome_melted$value,2)                    
ggbarplot(multiome_melted, 
          x = "library_name", 
          y = "value",
          fill = "variable", palette = "Paired",  label = TRUE,
  ggtheme = theme_pubr(x.text.angle = 45),
  position = position_dodge(0.9))

The fraction of sequenced reads that map to a unique gene in the transcriptome is expected to be higher than 50% (blue dashed line). However, the percentatge of reads mapped confidently to intergenic regions and reads mapped antisense to gene are expected to be lower than 30% (red dashed line).

mapping_qc_vars <- c(
  "library_name",
  "GEX.Reads.mapped.confidently.to.intergenic.regions",
  "GEX.Reads.mapped.confidently.to.intronic.regions",
  "GEX.Reads.mapped.confidently.to.exonic.regions",
  "GEX.Reads.mapped.confidently.to.transcriptome",
  "GEX.Reads.mapped.antisense.to.gene")

mapping_qc_gg.melt <- melt(qc_samples[,mapping_qc_vars])

ggplot(mapping_qc_gg.melt, aes(variable, value)) +
  geom_boxplot(fill = NA, outlier.shape = NA) +
  geom_jitter(width = 0.1, height = 0) +
  theme_minimal() +
  theme(legend.position = "none",
      axis.title = element_text(size = 12),
      axis.text = element_text(size = 12, color = "black")) + coord_flip() + 
  geom_hline(yintercept = c(0.3,0.5), col = c("red","blue"),  linetype='dotted')

6 Sequencing Saturation

Thirdly, we will plot the number of detected genes per library as a function of the total reads sequenced. We know that this function reaches a plateau, whereby more sequenced reads does not result in more detected genes. In those scenarios, we say that we have sequenced until saturation:

ggplot(qc_samples, aes(GEX.Sequenced.read.pairs, GEX.Total.genes.detected)) +
    geom_point() +
    geom_text_repel(data = qc_samples, aes(label = library_name), color = "black") +
    labs(x = "Number of Reads", y = "Total Genes Detected", color = "") +theme_classic() +
    theme(axis.title = element_text(size = 13, color = "black"),
          axis.text = element_text(size = 12, color = "black"),
          legend.text = element_text(size = 12, color = "black"))

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS/LAPACK: /Users/pauli/opt/anaconda3/envs/Tonsil_atlas/lib/libopenblasp-r0.3.10.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] plotly_4.9.2.1    ggrepel_0.8.2     kableExtra_1.3.1  knitr_1.30        data.table_1.13.2 reshape2_1.4.4    plyr_1.8.6        ggpubr_0.4.0      ggplot2_3.3.2     BiocStyle_2.16.1 
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5          lattice_0.20-41     here_1.0.1          tidyr_1.1.2         rprojroot_2.0.2     digest_0.6.27       R6_2.5.0            cellranger_1.1.0    backports_1.2.0     evaluate_0.14       httr_1.4.2          highr_0.8           pillar_1.4.6        rlang_0.4.8         lazyeval_0.2.2      curl_4.3            readxl_1.3.1        rstudioapi_0.12     car_3.0-10          Matrix_1.2-18       rmarkdown_2.5       splines_4.0.3       labeling_0.4.2      webshot_0.5.2       stringr_1.4.0       foreign_0.8-80      htmlwidgets_1.5.2   munsell_0.5.0       broom_0.7.2         compiler_4.0.3      xfun_0.18           pkgconfig_2.0.3     mgcv_1.8-33         htmltools_0.5.0     tidyselect_1.1.0    tibble_3.0.4        bookdown_0.21       rio_0.5.16          viridisLite_0.3.0   crayon_1.3.4        dplyr_1.0.2         withr_2.3.0         grid_4.0.3          nlme_3.1-150        jsonlite_1.7.1      gtable_0.3.0        lifecycle_0.2.0     magrittr_1.5        scales_1.1.1        zip_2.1.1           stringi_1.5.3       carData_3.0-4       farver_2.0.3        ggsignif_0.6.0      xml2_1.3.2          ellipsis_0.3.1      generics_0.1.0      vctrs_0.3.4         openxlsx_4.2.3     
## [60] RColorBrewer_1.1-2  tools_4.0.3         forcats_0.5.0       glue_1.4.2          purrr_0.3.4         hms_0.5.3           abind_1.4-5         yaml_2.2.1          colorspace_2.0-0    BiocManager_1.30.10 rstatix_0.6.0       rvest_0.3.6         haven_2.3.1